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We demonstrate that inverse statistical mechanical opti- 
mization can be used to discover simple (e.g., short-range, 
isotropic, and convex-repulsive) pairwise interparticle po- 
tentials with three-dimensional diamond or simple cubic 
lattice ground states over a wide range of densities. 

The properties of condensed phases are often linked to their 
structure. For example, heterogeneous materials with three- 
dimensional (3D) dielectric diamond morphologies can ex- 
hibit a photonic band gap ' , making them useful architectures 
for applications that range from lasers and sensors to solar 
cells. Although alternative methods for fabricating such ma- 
terials have been recently introduced, considerable interest 
remains in understanding how to create systems that spon- 
taneously self-assemble into structures with desirable prop- 
erties. Moreover, since various aspects of the effective in- 
teractions between nanometer- to micron-scale particles can 
be tuned experimentally via modification of solution or parti- 
cle properties''^, the following fundamental materials design 
question becomes especially relevant. Which types of inter- 
particle potentials provide a thermodynamic driving force for 
the particles to self-assemble into a given target lattice? 

Results from statistical mechanical theories, computer sim- 
ulations, and experiments have produced valuable insights into 
how to design interparticle interactions for self-assembly into 
periodic structures. For example, it is widely appreciated that 
spherical particles with steeply repulsive interactions spon- 
taneously assemble into highly-coordinated 3D structures''^, 
such as the face-centered cubic (FCC) lattice, at sufficiently 
high particle concentrations. Interactions that favor a tar- 
geted low-coordinated lattice ground state over other com- 
peting structures can also be designed by introducing specific 
types of complexity into the interparticle potential (e.g., multi- 



wide range of densities-remains an interesting open question. 

Inverse statistical mechanical methods such as those pi- 
oneered in recent years by Torquato, Stillinger, and oth- 
can be used to address this question. In fact, focus- 
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pie wells ' , non-spherical particle shapes or orientation- 
dependent "patches" on a particle surface ^'~^^), but those 
phases are generally stable over narrow ranges of thermody- 
namic conditions '^'^^. On the other hand, whether interac- 
tions with considerably simpler functional forms can also pro- 
duce targeted low-coordinated 3D ground states-stable over a 
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ing on the specific case of two-dimensional systems, Marcotte, 
Stillinger, and Torquato ^^'^' have employed inverse design 
principles to discover isotropic, convex-repulsive potentials 
with low-coordinated square and honeycomb lattice ground 
states. In the present study, we build upon that insightful body 
of work to search for simple pair potentials with specific low- 
coordinated 3D lattice ground states that are stable over a wide 
range of density. We choose the symmetric Bravais simple cu- 
bic lattice and the asymmetric non-Bravais diamond lattice as 
our target structures. 

The pair potentials we consider in our optimization are 
isotropic, convex-repulsive, twice continuously differentiable, 
and short-ranged. They are described (in terms of a charac- 
teristic energy scale e and length scale a) by the functional 
form,y(r/o) =e{/(r/o)+/,hift(r/o)}//[(rcut-r)/o]. Here, 
/(r/a)-motivated by a recently introduced model -is given 
by 

H is the Heaviside step function, and /shift(''/<^) =X{r/G)^ + 
Y{r/a) +Z. The constants X, Y, and Z are impUcit func- 
tions of the other parameters in the potential via the con- 
straints, y(rcut/a) = y'(rcut/a) = y"(rcut/o) = 0. In this 
study, we set rcut/<J = 2.25. As a result, y(r/o)/e depends 
on eight dimensionless parameters (A, n, A-i, k\, 5], X2, ^2, 
82), but one of these is not free because we further require 
/(l) + /shift(l) = V(l)/e = l. 

We obtain optimized potential parameters for specific tar- 
get structures using a standard simulated annealing algorithm 
(e.g., as described in Corana et al. ^^). Our optimization goal is 
to maximize the range of density over which the target lattice 
is the ground state for the potential. One practical way of ac- 
complishing this is to maximize the number n of uniformly 
spaced densities within a wide range [pi,pi + (m — l)Ap], 
where < n < m, for which the zero-temperature chemical 
potential (molar enthlapy) of the target structure is lower than 
those of the competing periodic structures at the correspond- 
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ing pressures. 

To get a more concrete sense of the optimization problem, 
consider the /* cycle at a given simulated annealing tempera- 
ture TsA- We first choose the eight primary potential parame- 
ters randomly, and subsequently solve for constants X, Y, and 
Z via the aforementioned constraints. If these 1 1 parameters 
are inconsistent with a convex repulsive potential, we reject 
them and try again; otherwise, we rescale the potential (and 
hence the parameters A, A-i, A-2, X, Y, and Z) by the constant 
factor required to ensure /(a) +/shift(o) = 1. For the resulting 
i* trial potential, the zero-temperature pressure and chemical 
potential of the target lattice at density pi are computed. This 
chemical potential is then compared to that of all other lat- 
tices in the competitive pool of structures (discussed below) at 
the same pressure. Similar comparisons are also carried out at 
pressures corresponding to the other m — 1 target lattice densi- 
ties in the range of interest [pi + Ap, pi + (m— l)Ap]. The 
number of state points in this set for which the target structure 
has the minimum chemical potential in the competitive pool 
is labeled The minimum chemical potential difference be- 
tween the target and its competing structures considering all m 
pressures is labeled A/j, (a quantity which is negative if the tar- 
get lattice is favored for at least one state point; i.e., if n, > 0). 
We define the simulated annealing energy for the trial po- 
tential Eg^ as 

E'sA = -«i + H(A^i) [«,■ + exp(A^,ye)] (2) 

In other words, the trial potential will be accepted as the 
cycle's pair potential with the standard Metropolis probabil- 
ity min [l , exp ( {E'^ — -Esa) / -^sa) ] where is the Boltz- 
mann constant. To fully explore the parameter space, we carry 
out optimizations initialized with various simulated annealing 
temperatures, and potential parameters. 

To ensure success of this optimization strategy, the com- 
petitive pool should ideally consist of all lattices which have 
chemical potentials that are similar to (or less than) that of 
the target structure for the class of pair potentials and state 
points under consideration. Motivated by the results of an ex- 
tensive ground-state study on related models we choose lat- 
tices for the competitive pool from the following types of pe- 
riodic structures: face-centered cubic (FCC), body-centered 
cubic (BCC), diamond (DIA), simple cubic (SC), wurtzite 
(WUR), hexagonal (SH), body centred orthorhombic (BCO), 
rhombohedral (hR), A7, A20 and pSn. Based on extensive 
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0.73 
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0.40 


5.32 


0.25 
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0.53 
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Table 1 Optimal parameters of the potential in equation 1 for 
diamond (DIA) and simple cubic (SC) target lattices. 
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Fig. 1 Ground-state phase diagrams of the optimized potentials 
(given in Table 1) for - (top) a diamond (DIA) target lattice, and 
(bottom) a simple cubic (SC) target lattice. Lattice parameters are 
reported in the supplementary material^ . 



preliminary calculations that we caiTied out for this study- 
which involved optimizing potential parameters using simu- 
lated annealing and computing ground state phase diagrams- 
we select the following specific lattices in the competitive 
pools (adopting previously introduced parameter nomencla- 
ture^') for use when the target lattice is diamond [FCC, WUR, 
SH (c/fl = 1.5), pSn (c/fl = 1.39), pSn (c/a = 1.25), A7 
= 3.79, M = 0.1385), A20 1.728, c/a = 0.626, 3; = 

0. 167)] and when the target lattice is simple cubic [FCC, BCC, 
DIA, SH (c/fl = 1), SH (c/fl = 1.08), SH {c/a = 1.172), A20 
{b/a = 1.72,c/fl = 0.66,3/ = 0.67), pSn {c/a = 0.873), pSn 
(c/fl = 0.78), pSn (c/fl = 1.75)]. Other lattices with differ- 
ent parameters may, of course, turn out to be more stable than 
these or the target structures under a given set of conditions. 
To test this possibility, we must carry out a more extensive 
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Fig. 2 (left) Optimized potentials (from Table 1 and equation 1) for diamond and simple cubic (SC) target ground states from this work; 
(center) optimized potential for the SC lattice from this work and a SC-forming potential developed by Fomin et al.^'^ ;( right) optimized 
potential for the diamond lattice from this work compared with an isotropic star polymer interaction model and the Yoshida-Kamakura 
potential which also exhibit diamond ground states. 



"forward" calculation of the ground state phase diagram with 
our final optimized potentials. 

There are a number of sophisticated search routines (e.g. 
genetic algorithms ■'^"■'^ and metadynamics-'^) that have been 
developed to find the most stable subset of lattices to con- 
sider in the forward calculation for a given potential. In this 
work, we construct the zero-temperature phase diagram for the 
optimized potential by searching for the most stable lattices 
from among the periodic structures mentioned above. For the 
structures which are defined by lattice parameters (WUR, SH, 
BCO, hR, A7, A20 and pSn), we use simulated annealing to 
obtain the optimal values of these parameters (i.e., those that 
minimize chemical potential) as a function of pressure. We 
also verify the mechanical stability of the optimized lattices 
on the phase diagram by analyzing their phonon spectra. The 
phase diagram is then constructed from among these energeti- 
cally and mechanically stable structures with optimized lattice 
parameters. 

Parameters of the pair potentials optimized for diamond and 
simple cubic target structures, respectively, in our simulations 
are reported in Table 1 . The corresponding ground state phase 
diagrams are shown in Figure 1. More information on the 
lattices is provided in the electronic supplementary material^. 
The first point to note is that both optimization strategies are 
successful in producing their target ground states over a wide 
density range. The stable density (po^) range for the diamond 
phase is 0.308 [1.213,1.521] and for simple cubic phase is 



0. 226 [1 .308, 1 .534] on the phase diagrams for their respective 
optimized potentials. These density ranges are considerably 
larger than those exhibited by the few other published models 
with isotropic potentials that can display these phases. 

Two relevant comparisons that can be made for the stabil- 
ity range of the diamond structure are to the coarse-grained 
center-of-mass star polymer interaction model developed by 
Watzlawek et al.^^ and to another model introduced by 
Yoshida and Kamakura^^'^". Although neither strictly satisfy 
all of the "simplicity" constraints of our optimized model po- 
tentials, they are simple no less and have been shown to exhibit 
stable diamond structures on their phase diagrams. Adopting 
the same non-dimensional representation of the present study, 

1. e., y(l) = e, the star polymer potential^^ (molecules with 
/ = 20 arms and e = ksT) has a stable diamond density range 
(po^) of 0.169, while the Yoshida-Kamakura^' potential has 
a diamond phase density range of 0.175. Both density ranges 
are roughly half of that exhibited by the potential optimized 
for the diamond structure in the present study. For the sim- 
ple cubic stiTicture, there are even fewer relevant comparisons. 
The one model that we are aware of exhibits a stable simple 
cubic ground state in a narrow density range of 0.03 (w 13% 
of the range displayed by the optimized potential presented in 
this work). For comparison, we plot the optimized potentials 
from our study along with the other potential models discussed 
above in Figure 2. 

Although the primary focus of the present work is design- 
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Fig. 3 Evolution of the configurational energy during a canonical 
MC simulation at r=0.040 and p=1.35. The red circles show the 
evolution of a system initiated as a perfect diamond lattice. The blue 
and green curves (represent a total of 16 configurations initiated 
with a high-temperature fluid configuration) correspond to 
configurations that crystallize with and without defects, respectively. 
The inset provides the pair correlation functions for some of the 
final structures. The solid red curve represents the system initialized 
from a perfect diamond lattice, whereas the blue and green symbols 
correspond to systems initialized from the high-temperature fluid. 
Blue and green symbols correspond to configurations that crystallize 
with and without defects, respectively. 



ing target ground states stable over wide density ranges, we 
have also completed some Monte Carlo (MC) simulations to 
probe the thermal stability of the diamond-forming system in- 
troduced here. We first completed a series of canonical MC 
simulations with = 250 at p = 1.35 to examine the melt- 
ing and freezing behavior. To estimate the melting point, we 
allowed a diamond lattice to relax at several temperatures sep- 
arated by Ar = 0.005 (in units of z/k-g), and found that the di- 
amond lattice melted at temperatures of T = 0.075 and above. 
To better understand the assembly process, we allowed a liq- 
uid, initially equilibrated at T = 1 .0, to relax at several tem- 
peratures, and found that the system assembles into a diamond 
crystal at temperatures of T = 0.045 and below. Figure 3 
provides data related to this assembly process. Specifically, 
we show the configurational energy as a function of MC step 
for a system equilibrated at T = 0.040. In this case, we find 
that each of the 16 configurations examined crystallize during 
the simulation, with 8 of the configurations forming a defect- 
free diamond lattice and 8 of the configurations assembling 
into defective diamond crystals. The nature of the underly- 
ing lattice was verified by examining the pair correlation func- 
tions. Similar results were obtained with a system consisting 



ofN= 1024 particles. Collectively, these data suggest that the 
diamond system exhibits a first-order melting transition at p = 
1.35. 

We are now employing free energy MC methods to con- 
struct phase diagrams for the systems introduced here. Fig- 
ure 4 provides initial data related to the diamond-fluid satu- 
ration curve. These points were located by finding the tem- 
perature at which the Gibbs free energy of the fluid matched 
that of the diamond crystal along a given isobar. The temper- 
ature dependence of the fluid's Gibbs free energy was com- 
puted via a combination of grand canonical transition matrix 
MC^^ and isothermal-isobaric temperature expanded ensem- 
ble MC^^ simulations. The temperature dependence of the 
crystal's Gibbs free energy was computed via a combination 
of Frenkel-Ladd MC^^ and isothermal-isobaric temperature 
expanded ensemble MC^^ simulations. Our results point to 
a concave-shaped diamond-fluid saturation curve within the 
temperature-pressure plane, with the temperature maximum 
located at approximately T — 0.065, where P = 10 and p = 
1.31. The heating and cooling simulations outlined above 
were completed at a density slightly beyond this maximum 
point. The Yoshida-Kamakura system exhibits a similarly- 
shaped diamond-fluid saturation curve with a lower maximum 
melting temperature^^ of T = 0.047. These results suggest 
that the thermal stability of the current model exceeds that of 
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Fig. 4 Phase diagram for the diamond-forming potential. The main 
panel provides the temperature-density plane. Circles and squares 
provide the saturated densities of the diamond and fluid phases, 
respectively. The inset provides the phase diagram within the 
temperature-pressure plane. Dashed lines are simply guides to the 
eye. The statistical uncertainty of the simulation data is smaller than 
the symbol size. 
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the Yoshida-Kamakura model. 

To summarize, our investigation shows that it is possible to 
use techniques of inverse statistical mechanical optimization 
to obtain simple pairwise interaction forms with targeted low- 
coordinated three-dimensional structures stable over a wide 
density range. We will be presenting a detailed study of the 
Monte Carlo free-energy simulation methods to compute the 
thermal stability (i.e., the temperature-dependent phase dia- 
grams) of the systems introduced here in a future publication. 
Additionally, we plan to study whether a systematic coarse- 
graining strategy (e.g., relative entropy maximization'*'') could 
preserve low-coordinated ground states when mapping from 
anisotropic "patchy" interactions to simpler, isotropic effec- 
tive potentials. 

As a final note, after finishing this manuscript, we became 

aware of a very recent preprint by Marcotte et al. ^' which also 
reports a convex-repulsive pair potential that exhibits a dia- 
mond ground state. They use a different inverse statistical me- 
chanical optimization method than that reported in this study 
with and obtain a considerably narrower range of thermody- 
namic stability. However, the resulting interaction potential 
and its derivatives, although of different functional forms, are 
strikingly similar to those we report here, providing further 
confirmation of the robustness of the qualitative result. 
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at Austin and the Center for Computational Research at the 
University at Buffalo for providing HPC resources that have 
contributed to the research results reported within this paper. 

References 

1 K. M. Ho, C. T. Chan and C. M. SoukouUs, Phys. Rev. Lett., 1990, 65, 
3152. 

2 A. Yethiraj and A. van Blaaderen, Nature, 2003, 421, 513. 

3 C. Likos, Phys. Rep., 2001, 348, 267. 

4 P. N. Pusey and W. van Megen, Nature, 1986, 320, 340. 

5 P. N. Pusey, W. van Megen, P. Bartlett, B. J. Ackerson, J. G. Rarity and 
S. M. Underwood, Phys. Rev. Lett., 1989, 63, 2753. 

6 M. C. Rechtsman, F. H. StiUinger and S. Torquato, Phys. Rev. E, 2007, 
75, 031403. 



7 M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. E, 2006, 
74, 021404. 

8 S. C. Glotzer, Nature, 2012, 481, 450. 

9 P. F. Damascene, M. Engel and S. C. Glotzer, Science, 2012, 337, 453. 

10 O. Gang and Y. Zhang, ACS Nana, 2011, 5, 8459. 

1 1 Zhang and S. C. Glotzer, Nano Lett., 2004, 4, 1407. 

12 E. Jankowski and S. C. Glotzer, Soft Matter, 2012, 8, 2852. 

13 H. Xiong, M. Y. Sfeir and O. Gang, Nano Lett., 2010, 10, 4456. 

14 F. Romano and F. Sciortino, Nature Comm., 2012, 3, 975. 

15 G. Doppelbauer, E. G. Noya, E. Bianchi and G. Kahl, Soft Matter, 2012, 
8, 7768. 

16 F. J. Martinez-Veracoechea, B. M. Mladek, A. V. Tkachenko and 
D. Frenkel, Phys. Rev. Lett., 2011, 107, 045902. 

17 E. G. Noya, C. Vega, J. P K. Doye and A. A. Louis, /. Chem. Phys., 2010, 
132, 234511. 

18 F. Romano, E. Sanz and F. Sciortino, J. Chem. Phys., 201 1, 134, 174502. 

19 S. Torquato, Soft Matter, 2009, 5, 1157. 

20 M. C. Rechtsman, F H. Stillinger and S. Torquato, Phys. Rev. Lett., 2005, 
95,228301. 

21 E. Marcotte, F. H. Stillinger and S. Torquato, Soft Matter, 201 1, 7, 2332. 

22 E. Marcotte, F. H. Stillinger and S. Torquato, J. Chem. Phys., 2011, 134, 
164105. 

23 H. Cohn and A. Kumar, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 9570. 

24 E. Edlund, O. Lindgren and M. N. Jacobi, Phys. Rev. Lett., 2011, 107, 
085503. 

25 Y. D. Fomin, E. N. Tsiok and V. N. Ryzhov, J. Chem. Phys., 2011, 134, 
044523. 

26 A. Corana, M. Marchesi, C. Martini and S. Ridella, Assoc. Comput. 
Mack, Trans. Math. Software, 1987, 13, 262. 

27 Y. D. Fomin, N. V. Gribova, V. N. Ryzhov, S. M. Stishov and D. Frenkel, 
/ Chem. Phys., 2008, 129, 064512. 

28 M. Watzlawek, C. N. Likos and H. Lowen, Phys. Rev. Lett., 1999, 82, 
5289. 

29 T. Yoshida and S. Kamakura, Prog. Theor Phys., 1972, 47, 1801. 

30 S. Kamakura and T. Yoshida, Prog. Theor Phys., 1972, 48, 21 10. 

31 S. Prestipino, F. Saija and G. Malescio, Soft Matter, 2009, 5, 2795. 

32 D. Gottwald, G. Kahl and C. N. Likos, / Chem. Phys., 2005, 122, 204503. 

33 T. Tuckmantel, R Lo Verso and C. N. Likos, Mol. Phys., 2009, 107, 523. 

34 E. Bianchi, G. Doppelbauer, L. Filion, M. Dijkstra and G. Kahl, J. Chem. 
Phys., 2012, 136, 214102. 

35 J. Behler, R. Martoiiffl£, D. Donadio and M. Parrinello, Phys. Rev. Lett., 
2008, 100, 185501. 

36 J. R. Errington, / Chem. Phys., 2003, 118, 9915-9925. 

37 A. P. Lyubartsev, A. A. Martsinovski, S. V. Shevkunov and P. N. 
Vorontsov-Velyaminov, y. Chem. Phys., 1992, 96, 1776-1783. 

38 D. Frenkel and A. J. C. Ladd, J. Chem. Phys., 1984, 81, 3188-3193. 

39 F. Saija, S. Prestipino and G. Malescio, Phys. Rev. E, 2009, 80, 031502. 

40 M. S. Shell, J. Chem. Phys., 2008, 129, 144108. 

41 E. Marcotte, F. H. Stillinger and S. Torquato, eprint arXiv:1212.3657vl - 
cond-mat.soft. 



5 



